A HAAR-TYPE APPROXIMATION AND A NEW NUMERICAL SCHEMA 
FOR THE KORTEWEG-DE VRIES EQUATION 

JASON BAGGETT, ODILE BASTILLE, AND ALEXEI RYBKIN 

I Abstract. We discuss a new numerical schema for solving the initial value prob- 

/*^ lem for the Korteweg-de Vries equation for large times. Our approach is based 

fvj upon the Inverse Scattering Transform that reduces the problem to calculating the 

reflection coefficient of the corresponding Schrodinger equation. Using a step-like 
approximation of the initial profile and a fragmentation principle for the scattering 
data, we obtain an explicit recursion formula for computing the reflection coeffi- 
cient, yielding a high resolution KdV solver. We also discuss some generalizations 
of this algorithm and how it might be improved by using Haar and other wavelets. 
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In this paper, we consider the well-known Korteweg-de Vries (KdV) equation 

Mf — 6uUx + Uxxx = 

on the whole line with the initial condition m(x, 0) = V{x). The function V{x) is 
assumed to be finite, nonpositive, and have compact support (i.e. zero outside of a 
finite interval). In particular, we will discuss a new algorithm for numerically ap- 
proximating the KdV for large times t. For small times, there are many available 

^jl algorithms to numerically integrate the KdV. Of particular interest is the opera- 

1/"^ tor splitting algorithm discussed in |[9|[l2). These two algorithms can be coupled 

C^ together to form a hybrid solver suitable for all times. 

r~^ The KdV equation is "exactly solvable" by relating it to the Schrodinger equa- 

O tion 

. . One obtains the so-called scattering data from the Schrodinger equation, and then 

the solution to the KdV can be obtained by performing the Inverse Scattering 
Transform (1ST). In this sense, the 1ST linearizes the KdV (as w^ell as some other 
^ nonlinear evolution PDEs). This provides us with an extremely powerful tool to 

Cd analyze its solutions. Unfortunately, numerical algorithms based upon the 1ST are 

much less impressive and have not so far shown a noticeable improvement over 
conventional methods. The real power of the 1ST is in capturing the large time 
asymptotic behavior of solutions to the KdV equation (e.g. solitons) which is of 
particular interest in applications. As presented in IIJ, an asymptotic formula for 
the solution of the KdV can be obtained from the 1ST Using this asymptotic for- 
mula, the large time solution of the KdV can be approximated from the scattering 
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data alone without the full manchineiy of the 1ST. This gives us a faster and more 
accurate method than other conventional methods for studying the large-time be- 
havior of solutions to the KdV. Moreover, although our method is a PDE solver, it 
does not use any standard numerical PDE techniques; we need only calculate the 
scattering data using root finders and some linear algebra as described below. 

The scattering data of the Schrodinger equation consists of the finitely many 
bound states — k^,, the corresponding (left) norming constants c^, and the (right) 
reflection coefficient R. The bound states are precisely the eigenvalues A of the 
Schrodinger equation that give square-integrable solutions (p. The left and right re- 
flection coefficients L and R, respectively, and the transmission coefficient T come 
from the asymptotic behavior of the left and right Jost solutions to the Schrodinger 
equation (pi and ^j., respectively, where for A = k^ 



4>i{^'k) - ^ ^,,^ .^^^. 



e'*^ + L{k)e-'^'' + o(l) x ^ -a 
T(A:)e'*^-^' x ^ oo. 



<Pr{^,^)-^^n...-,k^ 



and 

'g-fa _^ K(A:)e*^ + o(l) x ^ oo 
T(A:)e-'*^ x ^ -oo. 

If A = — K^ is a bound state, then (p^ and (p^ are square-integrable. The correspond- 
ing left and right norming constants are defined by 

^l ~ [ j \<Pi{x,iK)T^^[iK)\^dx\ c,. = ( / |^|,(x, fK)T~^(fK)|^dx I 

With our assumption that V has compact support, the bound states can be ob- 
tained as poles of R in the upper-half complex plane, and the (left) norming con- 
stants can be retrieved from the residues at these poles. The poles of R can be 
numerically approximated by using root finders. However, computing residues is 
numerically more difficult. We will instead consider a related function B w^hich is 
a rotation of the left reflection coefficient L. Then B has the same poles as R, and 
its corresponding residues are equal to the residues of R times a computable scal- 
ing factor. We give a new algorithm for computing the residues of B as presented 
below. 

We will approximate our potential using N piecewise-constant functions. Then 
in each interval where the function is constant, the reflection and transmission 
coefficients L„, Rn, and T,, can be explicitly derived. Let 

/1/T -R/T 
\L/T 1/T 

Then the reflection and transmission coefficients L, R, and T for the total potential 
can be derived from the principle of potential fragmentation (see, e.g. ||2i|8J), or 
layer stripping as it is known in the context of the Helmholtz equation |21[ : 

A = AN...A2A1 

for A: £ ]R where bars denote complex conjugation and A„ are the transition matri- 
ces 

1/T„ -Rn/Tn 



Lyi/ in t/ Jn 
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This gives us a recursive formula for the left and right reflection coefficients and 
also for the function B. Using this recursive formula for B, we can derive a recur- 
sive matrix formula for the residues of B at the poles in the upper-half complex 
plane. Consequently, all scattering data can be obtained, and the solution to the 
KdV can be numerically approximated for large times by the asymptotic formula 
given in |T|. 

In this paper, we also provide some numerical simulations. In particular, we 
give a comparison of computing bound states as poles of R and B and computing 
norming constants with our algorithm as opposed to other common algorithms. 
Although our algorithm is slower than standard methods for obtaining the scatter- 
ing data, we demonstrate that it tends to be more accurate, especially for discontin- 
uous initial profiles. We do not provide any error estimates; instead, the accuracy 
is verified on explicitly solvable examples. We also provide a comparison of the 
asymptotic solution to the KdV versus numerically integrated solutions. 

Lastly, we also give some generalizations of our algorithm. For example, there 
is a natural generalization of our algorithm to higher order spline interpolants of 
V{x). We also discuss possible improvements by using Haar and other wavelets. 
The Haar wavelets are piecewise constant functions that form an orthogonal sys- 
tem. Wavelets are closely related to Fourier series, and they exhibit many prop- 
erties that are numerically desirable. Since we are approximating our potentials 
V{x) using piecewise constant functions, one would believe that our algorithm 
can be modified to use Haar wavelets (and possibly more general wavelets). 

2. Notation 

We will denote the upper-half complex plane by C+. For a function / : C ^- C, 
we will let /(z) denote complex conjugation and f{z) = /(— z). As is customary 
in analysis, we will let L^{^) be the class of functions / such that L |/p < oo. 
We will let L{(]R) denote the class of functions/ such that /jj^(l + |x|)|/(x)| < oo. 
Given functions f,g& L^ (K), vve define the L^ inner product to be 



if-sh-Jj-s. 



and the L^-norm [| ■ [| 2 to be the norm with respect to this inner product, i.e. 

Il/ll2= / I/I' 

For a set A C C, Xa ^^i^ denote the characteristic function on A; i.e. 

f 1 iix e A 
I otherwise. 

3. Direct/Inverse Scattering Theory for the Schrodinger Equation 

ON THE Line 

Consider the KdV equation 

Ut — buUx + Uxxx ~ 0. 
A particular stable solution of the KdV equation is given by 

u{x,t) = -Ik'- sech^ (kx - ^}C't + 7) 
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where k and 7 are real constants. Solutions of this form are called solitons. For 
more general initial profiles u(x,0), in order to solve the KdV equation, we must 
first consider the Modified KdV equation (mKdV) 

Vt - 6v^Vx + Vxxx = 0. 

Miura (1967) discovered that one could obtain a solution to the KdV from a solu- 
tion of the mKdV by the transformation u = Vx + v^. Using Miura's transforma- 
tion, we obtain 

Ut - 6uUx + Uxxx = hv + — \ {Vt - 6v^Vx + Vxxx)- 

Through translation and scaling, one can transform the equation 

Ut — 6uux + Uxxx + Am = 
into the KdV 

Ut - 6uUx + Uxxx = 0. 
Because of this, Miura's transformation takes the more general form m — A = 
Vx + v'^. If we assume that ^ = %, then Miura's transformation u — \ = Vx + v^ 
becomes 

-^xx + ^i^'i)^ = ^^- 

This equation is fundamental to Quantum Mechanics. It is known as the one- 
dimensional, time-independent Schrodinger equation, where (p is the wave func- 
tion, M is the potential, and A is the energy. The problem of solving the KdV equa- 
tion reduces to finding nontrivial solutions of the Schrodinger equation in L^(]R). 
However, not every A has such a solution. Hence, the Schrodinger equation is 
an eigenvalue problem. Moreover, the eigenvalues of the Schrodinger equation 
do not change over time. For this reason, we can replace u{x,t) with our initial 
profile m(x,0) = V{x), and solve 

-^xx + V{x)cl> = Acp 

Suppose that V £ L| (M) . This ensures that there are finitely many soliton solu- 
tions. We have then that V{x) ^^ as x ^^ ±00. Hence, our solutions behave 
asymptotically like 

-'Pxx ~ H 

Since cp is bounded, we must have that cp behaves asymptotically like a sinusoid 
(A > 0) or decays like an exponential function (A < 0). 

Let H = —j-2 + V{x). Then the Schrodinger equation becomes Hep = Acp. The 
operator H is called the Schrodinger operator. We have that A is an eigenvalue of 
the Schrodinger operator H if H — A has no bounded inverse, and we say that A 
is in the spectrum of H. For each A > 0, there is a nontrivial solution to Hep = Acp. 
However, these eigenf unctions cp are not contained in L^(]R). We call the set of 
such A the continuous spectrum of H. The eigenvalues A < give square-integrable 
eigenfunctions cp. However, there are only finitely many such A. We call these A 
the bound states, while the set of bound states is called the discrete spectrum. The 
continuous spectrum gives rise to a component of the solution of the KdV which 
acts like a solution to the linear equation Ut + Uxxx = 0. This part of the solution is 
the dispersive portion of the wave. The discrete spectrum corresponds bijectively 
with the soliton solutions. This portion of the solution of the KdV is stable and 
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(a) Wave e^'^^ radiating from oo 




(b) R(*:)e*-"' is reflected 
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(C) T{k)e^'" is transmitted 




T(k)e" 



(d) Similarly, we can consider a 
wave coming from — oo 



Figure 1. The Jost solutions (pi{x) and ^,.(x) as waves radiating 
from ±00 



does not decay over time. Thus, we are really only interested in knowing the 
discrete spectrum for large times. 

Suppose A = A:^ G R. Among all solutions (p to the Schrodinger equation, pick 
the one satisfying 



(p,{x,k) = 



-ikx 



R(A:)e"^^ + o(l) x ^ oo 



T{k)e- 



-ikx 



X — > —00. 



Such a solution (p^. is known as the right Jost solution. The function T{k) is known as 
the transmission coefficient, and R{k) is the right reflection coefficient. The reason for 
this terminology is that we can view (p^. as a wave e~ radiating from infinity, and 
R{k)e is the portion of the wave that is reflected while T{k)e^ is the portion 
that is transmitted (see FigurellJ. Similarly, we can consider the left Jost solution 



4>i{x,k) 



e'^'' + L{k)e-'^^ + o{l) x ^ -oc 
T(A:)e'*^-^ x ^ oo. 
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where T{k) is the same transmission coefficient and L{k) is the left reflection coeffi- 
cient. 

li k = ix is a bound state, then ^i{x,iK),(pj.{x,iK) e L-^(]R). We define the left 
and right norming constants atk = ix to be 

C; = \\(pi{x,iK)T^'^{iK)\\2^ Cr = H^^ (x, /k)T"^ (k) Hj"^ 

4. The Classical Inverse Scattering Transform 

Since V{x) G LJ(]R), we have that there are finitely many bound states A = k^ 
where k = ix. Let K denote the number of bound states, and let 

Ki > K2 > ... > Kk> 0. 

Let c„ denote the left norming constant at A: = iK„ . 

Once we know the scattering data for the Schrodinger operator, we can use the 
Inverse Scattering Transform (1ST) to obtain the soliton solutions of the KdV equa- 

.• , , direct scatterine , , 

tion. u{x,0) ^^ S(0) 

time evolution 

u(x,t) ^ Sit) 

inverse scattering 

The Direct Scattering Transform is we map the initial potential m ( j, 0) into the scat- 
tering data 

Next, we evolve the scattering data over time in a simple fashion: 

• C„{t) = Cne^^nt 

• R{k,t) =K(fc)e8*'f 

where k„ = Kn(0), c„ = c„(0), and R{k) = R{k,0). Then the scattering data 
becomes 

S{t) = {{-K„{tf,c„{t)}^^^,R{k,t),k e R}. 
We can reclaim the solution to the KdV using Inverse Scattering as follows: 

• Form the Gelfand-Levitan-Marchenko (GLM) kernel: 

F(x, = E cl{t)e-''"'' + -^ / e*^R(fc, t)dk. 

• Solve the Gelfand-Levitan-Marchenko equation for K{x,y, t), y > x: 

/"CO 

K{x,y,t)+F{x + y,t)+ / F(s + y, t)K{x, s, t)ds = 0. 

• The solution to the KdV equation is 

u(x,t) = —2z—K(x,x,t). 
dx 

Luckily, for large times t we can simplify the GLM kernel. We have that 

e'''''R{k, t)dk ^ as t ^ oo. 
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Thus, for large times we can approximate the GLM kernel by 



N 






Let 



C(x,0) = 


/Cii(x) Ci2(x) 
C2l{x) C22{x) 


... Ci^n(^)\ 
... C2,n{x) 


\Cn,i(x) Cn^2{x) 


■■■ cn,n(x)/ 


where 




Cmn[X) = e ^ 

Km + Kfi 


Kui-t-KjijX 


The matrix C evolves in time 


by 





Then for large times, our solution to the KdV is 1 1 1 

u{x,t) w -2-^ln[det(I + C(x,0)] 
From this, one obtains the asymptotic formula |jlj 

N 



u{x,t) -^ — 2 ^ K^, sech2(K„x — 4K^f + ^1-^7^) 



where 



Notice that the large time solution u{x, t) of the KdV behaves like a finite sum of 
single solitons. Moreover, we no longer need to do the full 1ST to solve the KdV 
for large times. We need only find the bound states — k^, and norming constants 

Cn- 

If R and T can be analytically continued, then the poles of R and T in C+ are 
precisely iK„. That is, all of the poles of R and T in C+ lie on the imaginary axis 
and correspond with the bound states. Better yet, these poles are actually simple 
p)[T5). Furthermore, if we assume that V has compact support, then ||5l 



Res R{k) = icl. 

k=iKn 



(4.1) 



Consequently, if our potential V has compact support and R is analytically contin- 
ued into C"'", then the bound states and norming constants can be obtained from 
knowledge oi R{k) for k e C^. In this case, we can approximate the solution of 
the KdV for large times from only knowledge of R{k) for k E C+. 

5. The scattering Quantities for a Block (Well) Potential 

Consider the case w^hen our potential V is a single nonpositive well w^hich is 
—a^ on the interval [—b,0] and elsewhere, i.e. V{x) = —a^x\-b,o]{^) (^^^ Figure 
I2I. In this case, we can obtain an exact solution to the Schrodinger equation. More- 
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Figure 2. The setup for a single block potential 



over, using the continuity of the solution ^ and its derivative (p^, we can set up a 
system of equations and solve for R and T. Doing this, we obtain 



R{k) = w 
where 






-iab{a>—l/a>) 



1 4 

^n \ -1- — <^ iit 



(5.1) 



a; 




c.iX'i'+l/a;)!;?; 



These formulas for R, L, and T are actually meromorphic in C if we choose the 
branch cut along the imaginary axis between —ia and ia. Using these formulas, R, 
L, and T can be analytically continued in C+. The only difficulty lies in considering 
the branch cut. However, we have that ci>{—k) = —co{k) and ^{—k) = ^{k)- It 
follows that R{-k) = R{k) and T{-k) = T{k). For k e flR, we have that R{k) = 
R{-k) = R{k) and T{k) = T{-k) = T{k), so R and T are real-valued for k e flR. 
For A: = +0 + zk, we have that —k= — + zk. Therefore, R{—0 + /k) = R(+0 + ix) 
and since R is real-valued on flR, R(+0 + ix) = R{+0 + ix). Hence, R(— + zk) = 
R(+0 + zk), so R is continuous along the branch cut between —ia and ia. It follows 
that R is meromorphic in C Similarly, T is meromorphic in C as well. 

Consider the poles iK„ and residues ic„ of R. The poles of R and T satisfy ^ = 
a;4. If we let Vn = ^, then k„ and c„ can be explicitly computed by the following 
formulas: 



ah 

n 



1 



n 



arctan 



1 






(5.2) 



and 



for n = 1,..., 



2k„ (l - (^)') 



2 + ;7K„ 



(5.3) 
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6. The Potential Fragmentation and the Scattering Quantities for 
Potentials Composed of Blocks 

We define the scattering matrix to be 

'T R 



The matrix S is unitary, i.e. S^^ = S* where S* is the conjugate transpose of S |3j. 
This gives us a few identities, namely |6{15J 

Lf + TR = (6.1) 

for A: £ ]R. If we were to shift our potential to the right by p, then the scattering 
matrix would change as follows Q: 

L{k -p)= L{k)e^'''P (6.2) 

T{k - p) = T{k) (6.3) 

R{k -p)= R[k)e-^^^ (6.4) 

Now suppose that our potential V consists of N nonpositive blocks. Let Rn,L„, Tn 
be the reflection and transmission coefficients on the n-th block: Vn{x) = —a^ on 
[—bn, —bn-i] where Bq = 0. Let R^, L^,, T° be the reflection and transmission coef- 
ficients on the n-th block shifted to the origin: V„{x) = —a\ on [—{hn — fc„_i),0]. 
Let Ri,2,...,nf ^i,i,...,n, T^i,i,...,n be the reflection and transmission coefficients on the 
first n blocks. R, L, T without subscripts or superscripts will denote the reflection 
and transmission coefficients for the overall potential. 

Let 

fl/T -R/T\ 

The fragmentation principle, or layer stripping principle as it is also known, p}[6p)[2T| 
says that for A: G ]R 

A = AN...A2A1 (6.6) 

where A„ are the transition matrices 

^ fl/Tn -Rn/Tn' 

" \un„ 1/T„ 

(Note that blocks with fl„ = may be simply ignored since this implies A„ is 
the identity matrix). We can use potential fragmentation to come up with some 
recursive formulas. Using 6.1 6.2 [6.3 6.4 and 6.6 we obtain that 



^1 "+1 - " ^^ i_kO u ^2^ ■ ^^ 

A similar expression may be obtained for the left reflection coefficient: 

dO I „2ikh„ _ ^O"' 

^l,...,n+l — --^— 1 oO T Jilch ' ^ ' 

n+1 

We have that Li,...^„ = - JieS-K^ for A: e R. Thus, |Li,...^„| = |Ki,...,„| for A £ R. 

Ti,...,™ 



Therefore, Li^...,„ = Ki,...^„e"2(7c/J„ f^j. g^j^g /3„ : IR ^ R. Equations l|67|| and l|6^ 
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then give us that 

Rl,...,n+1 

where j6^ = bi and 



Ki 



^0 2/*:{fc„-/3J _ o 



^n+1 



K+lKZr^^ „.2,7:(fc„-/3„) ~ 






Ki 



R 



Define A„ 



( |6.10| give us that 



and 



i^u^^iikbn^ Then A„ = e2*(fe„-/3„) f^j. ;^ 



:|V+ig-2ifcfc„ 



(6.9) 



(6.10) 



G R. Equations | |6.9) and 



R 



l,...,n+l 






^nK+1 



Ri 



A 



n+l 



;i-A.K°+i^i n 

K+lKZr^nRl n-K^^ 






(6.11) 



(6.12) 



where Ai = 1. Let us next define B„ 
following recursive formula: 

Bn+l = 



^1 " A„RO^j - Ki „ 

AnRi,...,n = U,...,ne^'^^"- Then we get the 









(6.13) 



where Bj = i^i. Notice that B„+i is a Mobius transform of B„, and that the recur- 
sive formula for B„ is much simpler than the recursive formula for i^i^.. „ . More- 
over, this formula only depends on B„, i^°_|_i, and i^°+i- From ! 



5.1 



R^ 



^„ 



■07^ 



(6.14) 



where ft,, = b„ — b„_i is the width for the n-th block, a;,, (A:) = ^ + y (^)^ + 1/ 
and ^„(A:) = £->««'>« {a'„(/:)+l/a;„(»:))^ poj. fc g ]R^ we have that K^" = ^. By taking 
the complex conjugate of 1 6.14 1, we obtain that for A: G R, 



R° 



^n'^n 



(6.15) 



Since Kf, is meromorphic in C, Rn is meromorphic in C as well (in particular, both 
are meromorphic in C+ w^here the poles of interest lie). Since the formula in | |6.15| 
is meromorphic in C, it follows that l|6.15| holds for all A: £ C Continuing induc- 
tively using equations ( |6.11| , ( 6.12) , | |6.13) , it follows that Ki ...,„, A„, and B„ can be 
continued to meromorphic functions in C (and in particular, C^) for all 1 < « < N^- 
Since B„ = A„_Ri „ = Li „e^ «, we have that B„ has the same poles k = iK^ 
in C+ as Lj „. Consequently, B„ and Ki ...,« have the same poles in C+. Since the 

poles k = ixm in C+ of Li,...,„ and -Ri,...,„ are simple, we have that A„ = .^^^i^uHe^''^''" 
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is analytic in C+ and nonzero at all k = ixm- It follows from equation ( |4.1) then 
that 

Res B„ = An{iK„) Res Ki „ (6.16) 

= icl,An{iK„,) (6.17) 



The value of A„{iK,n) can be determined via the recursive formula 6.12 If we 
can determine an algorithm for determining the residues of Bjv, then this would 
effectively give us an algorithm for calculating the (left) norming constants. 
Now suppose B„ has the form 

B„ = -S^. (6.18) 

Applying | |6.13| , we get a linear system of recurrence relations for p„ and qn'- 

(6.19) 



Pn+l = -RnPn " R°+iRn'?n 

qn+1 = K+lKPn + R^nqn 



or in matrix form 



where 



P"+M = Mn (l") = M„...M2Mi ( PM (6.20) 



Let N denote the number of blocks. If qN{k) = but k is not a pole of Bjv, then 
Pn = as well. From ( 6.20| , this means that det(M„) = for some 1 < n < N — 1 



or I I = 0. Since Bi = Ri, from 1 6. 18) we have that j- = —Ri- Our choice of 



qij ^ - ' ' <?! 



pi and qi is arbitrary, so long is this ratio is preserved, since our resulting solution 
of B„^i is independent of our choice for pi and qi. Some choices for our initial 
vector may be preferable for numerical computations, but for our purposes we 

Pl\ _ (-R 



will choose ( I = ( i ^ ) ' because it is nonzero for all k. Thus, if ^n =0 but 

k is not a pole of Bjv, then det(MN-i...M2Mi) = 0. Equivalently if (Jn = and 
det(Mjv-i-M2Mi) 7^ 0, then A: is a pole of -Ri,...,n- 

We claim that det(MN-i---M2Mi) (A:) = for some A: e C+ if and only if A = mn 

or for some 1 < n < N — 1 and some Q < m < 



aJi 



. We have that det(Mjv-i...M2Mi) = if and only if det(M„) = for some 
1 < n < N — 1. Moreover, 



det(M„) = RORO(R°+iRO+i-l). 
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Thus, det(M„) = if and only if RO = (equivalently K° = 0) or K°_^iK°_^i = 1. 
The second case occurs when 

After some algebra and noting that ^„j^i = e~"^n+\K+i(^n+i+'^/'^n+i) ^ g, this sim- 
plifies to (^^ 1 2 = 1. A simple calculation then gives us that (-O^^^^ = 1 for A: G C"*" 



if and only if A: = Zfl„+i. After a lengthy computation using 1 6.14 1, we obtain that 
R° {k) = for A: e C+ if and only if equation | |6!2l) holds. 

Now suppose that ^^(fc) = 0, det(Mjv-i...M2Mi)(A:) 7^ at A: = z'k, and that k 
is not a pole of R^. Then A is a pole of Bjv, and consequently a pole of R = Ki^...,]v 
as well. Consequently, A'^ is a bound state of the Schrodinger equation. Since 
det(MN_i...M2Mi)(A) 7^ 0, we have that piv(^) 7^ 0- Since A is not a pole of R^ 
and since R^ and R^ have the same zeros with the same multiplicity, — -^p,, 7^ 0. 
However, q^^ = and the poles of Bjv are simple, so 

ResBN = -S^. (6.22) 



To find (j^, we can differentiate 16.19 1 to acquire 



^"+M=M„KM+M;K" . (6.23) 

Hn+l/ VI n/ \HnJ 

Therefore, for the poles of R w^here det(MN...M2Mi) 7^ and that are not poles 
of R^, the residues of R can be recovered through | |6.22) and | 6.17| . 



7. Numerical Simulations 

Tables [TII3I and [5] give a comparison of some algorithms for calculating bound 
states described below. The exact bound states in table IT] were calculated using 
equations | |5.2) . All calculations were performed using MATLAB. 

There are two commonly used numerical methods for approximating the bound 
states: 

(1) Matrix methods - Estimate the Schrodinger operator H = — jr + ^(^) 
using a finite-dimensional matrix and find the eigenvalues of the matrix. 
In particular, |22| describes how this can be done using the Fourier basis. 
In tables[l]|3] and|5] a 512 x 512 matrix is used. 

(2) Shooting Method - The Shooting Method involves recursively choosing 
values of A and "shooting" from both end points to a point c £ [a, h] . De- 
fine the miss-distance function to be the Wronskian determinant 



D(A) 



ml(c,A) wr(c,A) 
Ml(c,A) 4(c,A) ■ 

where ml is the interpolated function from the left endpornt and mr is from 
the right endpornt. If A is an eigenvalue that satisfies the boundary value 
problem, then D(A) =0. For more details, see for example flTi. 

There are of course other existing methods for approximating bound states; see for 
example [10 11 13 Wj. However, we will only focus on these two. 
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If one approximates the potential using finitely many blocks, then we can use 
the following algorithms for estimating bound states: 

(3) Use the recursive formulas 16.11 1 and 1 6. 12) to find the bound states as 
zeros of 1/-Ri,...,n- 



(4) Similarly, one can use the recursive formula 16.13 1 to find the bound states 
as zeros of l/Bj^j. 

(5) Using | |6.20) , the bound states can be found as zeros of q^^. One must also 
check the values of k listed in | |6.21) where det(Mjv-i...Mi) = 0. 



Algorithm (1) seems to be the fastest of these algorithms, followed closely by (2). 
Moreover, algorithm (1) has great accuracy w^hen the initial potential is smooth. 
However, for discontinuous potentials, the Gibbs phenomenon severely hinders 
the accuracy of the algorithm. Moreover, the domain chosen seems to effect al- 
gorithms (1) and (2) greatly. On the other hand, tables [3] and [5] demonstrate that 
algorithms (3)-(5) are more robust when choosing the domain, with a smaller do- 
main being preferable for the amount of time. All of algorithms (2)-(5) rely on 
finding roots of some function, so inheritently all of these functions have all of the 
problems that root finders tend to have. For example, given a good initial approx- 
imation of a bound state, the root finder might diverge or converge to a different 
bound state. Furthermore, the bound states are known to cluster towards 0, which 
makes it increasingly difficult to accurately determine all of the bound states as the 
number of bound states increases. However, w^hen the root finders do converge, 
algorithms (3)-(5) are extremely accurate. Algorithms (3)-(5) also seem to be much 
slower than algorithms (1) and (2), with (5) being the slowest. 

In summary, the commonly used algorithms (1) and (2) for calculating bound 
states are much faster than the other algorithms. Moreover, algorithm (1) tends 
to be extremely accurate, especially when the potential is smooth. However, al- 
though algorithms (3)-(5) are much slower, they also tend to be very accurate, 
especially with discontinuous potentials. Moreover, these algorithms seem to be 
more robust when choosing the domain of the potential. 

Supposing the bound states have been calculated, tables [2] and |4] give a com- 
parison of some of the various algorithms for computing (left) norming constants. 
First is the algorithm described in the present paper: 

(i) The potential is approximated using finitely many blocks, and the norming 
constants are calculated as residues via equations ( 6.22| and 1 6. 17) . 



Next we have the obvious algorithm using the definition of the left norming 
constant: 

(ii) Suppose V has compact support [A, B] . We have that (p{x, k) = cpf (x, k) /T{k) 
satisfies <p{x,k) = e for x > B. One can numerically integrate the 
Schrodinger equation from B to A. Then c^ = \\(p\\2 > ^vhich can be nu- 
merically integrated. 

The authors were also presented the following algorithms by Paul Sacks: letting 
fl = l/T and h = —R/T, then R = —^ and the transition matrix A given in | |6.5| 
becomes 
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Table 1. V{x) = 
step size h = 0.01 



-4:X\-4fl] {x), domain chosen [—10, 10], spacial 



Algorithm 



Kl 



K2 



«3 



Relative Error 



Time (sec) 



Exact 
(1) 
(2) 
(3) 
(4) 
(5) 



1.899448036751944 
1.898826427139628 
1.899418261950639 
1.899448036751942 
1.899448036751949 
1.899448036751942 



1.571342556813314 
1.568514453040000 
1.572105829640451 
1.571342556813313 
1.571342556813312 
1.571342556813315 



0.876610362727433 
0.867505110670815 
0.872097420881459 
0.876610362727439 
0.876610362727428 
0.876610362727434 




0.003651829842877 
0.001749410414267 
0.000000000000003 
0.000000000000003 
0.000000000000001 



0.004355000000000 
0.126239000000000 
0.505034000000000 
4.168762000000000 
5.425778000000000 
10.268152000000001 



Moreover, b is analytic everywhere in C^, and the simple poles of T in C^ 
simple zeros of a. Consequently, | |4.1[ | gives us that 

2 -biiKn) 



are 



a' iK„ 



The derivative a' with respect to k can be approximated using the central difference 

a{k + r]/2)-a{k-r]/2) 



a'{k) 



V 



The question then becomes how one evaluates a(k) and b{k). Here are two ap- 
proaches: 

(iii) The potential is approximated using a finite number of blocks, and a and 
b are calculated using potential fragmention | |6.6| . The transition matrices 
are evaluated using equation | |5.1[ |. 
(iv) Supposing the potential has compact support [a, /5] , the Schrodinger equa- 
tion can be numerically integrated from a to fi with the initial conditions 
(p{a,k) = e'''"',(p'{a,k) = -ike-'^". Then cp{x,k) = cp^{x,k)/T{k), so for 
x> p 

(p{x,k) = a{k)e-'^'' -b{k)e^\ 
Consequently, a and b can be retrieved from 



(a{k)\ 




(p{ci,k) 
(p'{a,k) 



Algorithms (ii) and (iv) seem to be the fastest of these four algorithms. How- 
ever, they are also the least accurate since they require integrating the Schrodinger 
equation which is extremely sensitive to errors. Algorithms (iii) generally takes 
about half as long as algorithm (i). Algorithm (i) seems to be the most accurate 
for discontinuous potentials, w^hile algorithm (iii) seems to be the most accurate 
for smooth potentials. Moreover, the accuracy of algorithms (i) and (iii) increases 
when the bound states are approximated using algorithms (3)-(5). 

Lastly, figures [3] and |4] compare the asymptotic formula given in fll with the 
numerically integrated solution obtained by using the split step Fourier method. 

8. Haar Systems and a KdV Large-Time Solver 

Suppose now that V is finite, nonpositive, and has compact support. Then V 
can be well approximated using finitely many nonpositive blocks. For such poten- 
tials V , we now summarize the algorithm for solving the KdV for large times: 
• Approximate the potential V{x) using N nonpositive blocks 
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Table 2. V{x) = — 4;(;r_4 gi (x), domain chosen [—4,0], spacial 
step size h = 0.01, energy step size i] = 0.001, exact bound states 
used 



Algorithm 



Relative Error 



Time (sec) 



Exact 0.038798932148319 0.145167980693995 

(i) 0.038798932148326 0.145167980694058 

(ii) 0.141300713908832 0.293968570328614 

(iii) 0.038798937542783 0.145168027811526 



0.005992000000000 
0.000000000002272 2.008827000000000 



0.257227284424067 

0.257227284424741 

0.444084025980906 0.872538777834092 0.032151000000000 

0.257226712349713 0.000001926938416 2.070128000000000 



(iv) 0.051311576782601 0.109225786002665 -0.041977058580690 1.012467614318974 0.147137000000000 



Tables. V{x) = — sech (x), domain chosen [—5,5], spacial step 
size h = 0.01 



Algorithm 


K 


Relative Error 


Time (sec) 


Exact 


1.000000000000000 








(1) 


1.000181385743159 


0.000181385743159 


0.123699000000000 


(2) 


1.000010661550817 


0.000010661550817 


0.165820000000000 


(3) 


0.999997769556372 


0.000002230443628 


8.624536000000001 


(4) 


0.999997769556371 


0.000002230443629 


8.780760000000001 


(5) 


0.999997769556372 


0.000002230443628 


14.264264000000001 



Table 4. V{x) = — sech^(x), domain chosen [—5,5], spacial step 
size h = 0.01, energy step size t] = 0.001, exact bound state used 



Algorithm 


c^ 


Relative Error 


Time (sec) 


Exact 


2.000000000000000 








(i) 


2.004086813877857 


0.002043406938928 


3.460512000000000 


(ii) 


1.274509474591987e-004 


0.999936274526270 


0.023956000000000 


(iii) 


1.999683571279579 


0.000158214360211 


1.631856000000000 


(iv) 


1.993946894122799 


0.003026552938601 


0.088449000000000 



Table 5. V{x) = 
step size h = 0.01 



sech (x), domain chosen [—10,10], spacial 



Algorithm 



Relative Error 



Time (sec) 



Exact 
(1) 
(2) 
(3) 
(4) 
(5) 



1.000000000000000 
1.000000008244449 
1.000071226867798 
0.999997777799307 
0.999997777799307 
0.999997777799307 




0.000000008244449 
0.000071226867798 
0.000002222200693 
0.000002222200693 
0.000002222200693 





0.119652000000000 

0.177487000000000 

17.526257000000001 

18.236691000000000 

28.798037999999998 



Bound states are found as zeros of l/i^i,...,]v with initial estimates, for ex- 
ample, derived from a spectral matrix estimate of the Schodinger operator 
The norming constants are calculated as residues of Bjv at the bound states 
using the previously described recursive formulas 
The solution to the KdV is obtained from the asymptotic formula: 



(x,i) 



N 



-2 J^ K^ sech (k„x ■ 

n=l 



■4K^f + ln. 



where 



7„ 



2Kn 



n 



There are a number of possible improvements to this algorithm. For example, 
the number of bound states is known for a single block, so the results in 17] could 
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Figure 3. V{x) = -lOsech^(x), t = 0.3 
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Figure 4. V{x) = -5sech2(x), t = 0.6 

possibly be implemented to determine the exact number of bound states for the 
potential. As another example, instead of piecewise-constant functions, one could 
instead use higher order spline interpolants of the potential. All of the recursive 
formulas in section [6] were derived from potential fragmentation, which holds for 
arbitrary potentials; the only things that would change would be the formula for 
R^, the initial values in the recursive formulas, and the values for k in (6.21[. For 
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example, in the case of piecewise-linear spline interpolants, the formula for R„ 
would involve the Airy functions. 

Another possible route for improvement would be the use of Haar wavelets or 
other wavelets. We will only consider Haar wavelets in the current paper. For a 



great exposition on Haar and other wavelets, see 1 16 1. Consider the scaling function 
(p{x) = (Pq{x) 



1 if < X < 1, 
otherwise. 



and the mother wavelet 

{1 if0<x<l/2, 

-1 ifl/2<x<l, 

otherwise. 

We form the Haar wavelets as follows: let 

Wjfi{x) = w(l]x). 

Then w-,q has support [0,2^^]. Next, we translate w^q so as to fill up the entire 
interval [0, 1] with 7) subintervals of length 2~^ 

Wj^l^{x) = cp2j^j^ = ^j,o{x ~^) ~ w{2^{x — k)), k = 0,1,. ..,2'' — 1. 

Then zfy jt has support [2~ik, 2~i (k + l)]. The collection of Haar wavelets 

^2« ={f„, ■■ 0<m<2"-l} 

forms an orthogonal system with respect to the L^ norm of dimension 2"; the col- 
lection "Hco forms a complete orthogonal system for L^([0, 1]). For 7i2", let (Pj. de- 
note the vector in ]R corresponding to (p/, i.e., the entries of q>y are the function 
values of (p^. on the 2" intervals. 

By translating and scaling, suppose without loss of generality that V has com- 
pact support [0,1]. Since V is finite, we have that V e L^([0, 1]), so V can be 
expressed in terms of the Haar basis: 

oo 
r=0 

where 



\\9r\\2 

Let Vq denote the piecewise-constant approximation of V on the 2" intervals 
mentioned above, and let V denote the corresponding column vector in IR^ . Then 
Vq can be represented as a linear combination of the Haar wavelets in 7^2" ■ 

2«_l 
r=0 

where the coefficients Cy are as described above. Letting c denote the column vec- 
tor of coefficients c,-, the discrete wavelet transform (DWT) is the map H2" : V l-^■ c; 
that is, H2" is a change of basis from the standard basis to the Haar basis. Letting 
W2« denote the matrix whose r-th column is cp^, we have that 

V = W2«C, 



18 UAF 



SO 



W2-„^V, 



implying that H21' = Wji,^ . (Note: often, the columns are normalized so that W2n 
is an orthogonal matrix. In this case, H2« = W|„ where * denotes the transpose). 

The Discrete Wavelet Transform is analogous to the Fast Fourier Transform 
(FFT), which expresses V in the orthogonal basis corresponding to the Fourier 
basis {e'^ ^ : —2"^^ < r < 2"} in L^([— tt, tt]). However, the Fourier basis is not 
localized, unlike the Haar basis, so the Fourier basis has difficulty capturing data 
concentrated in a relatively small region. The Fourier basis tends to accurately ap- 
proximate smoother functions, w^hile exhibiting the so called Gibb's phenomenon at 
discontinuities. On the other hand, the Haar basis tends to accurately approximate 
discontinuous functions, while only slowly converging to smoother functions. 

In the context of solving the KdV, Haar wavelets may possibly be implemented 
in a couple ways. One approach would be to approximate the potential using Haar 
wavelets since it generally gives more accurate piecewise-constant interpolants 
than, say the midpoint rule. Then the interpolating potential would be changed to 
the standard basis and used in our algorithm. 

Currently, our potentials are being approximated by step functions using the 
standard basis since this is the form required for potential fragmentation. How- 
ever, it is more desirable to represent the potential using Haar wavelets in many 
cases, such as for signal processing. Another approach for improving the algo- 
rithm would be to change all of our formulas over to the Haar basis. There are still 
many open questions in this regard, for example 

(I) If our potential was to be expressed in the Haar basis instead of the stan- 
dard basis, what would be an efficient way to determine the scattering 
data? 
(II) Could potential fragmentation and our recursive formulas be modified to 
use the Haar representation of the potential? 
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